Project: Defining new types



In [1]:

    
include("dual.jl")









    Out[1]:





log (generic function with 20 methods)

Let's implement a very basic version of dual numbers (done much better here. A dual number is written as $$a + b\epsilon$$ where $a$ and $b$ are real numbers and $$\epsilon^2 = 0$$ Basically, a dual number encodes both a point in space and a local derivative.

To start, we'll assume that dual numbers are only a generalization of Float64.

We will also define (though it's not rigorously true) that these numbers are subtypes of Real. That way, any function that admits a Real will also admit a Dual.

immutable Dual <: Real val::Float64 adj::Float64 end

Notes:

Again, we choose Dual <: Real for convenience, not math. We want Duals to go where Reals are allowed.
Instead of immutable we could have used type. immutable means that we cannot change these numbers, but as a consequence, they can be stack-allocated and faster.
We need to give Julia the types of the fields. (If we don't, they're assumed to be Any.) In contrast to function arguments, we absolutely want these to be concrete, if possible, since abstract fields lead to a performance hit.



In [2]:

    
aa = Dual(1, 2)  # note implicit conversion to Float64









    Out[2]:





1.0 + 2.0ϵ



In [3]:

    
fieldnames(aa)









    Out[3]:





2-element Array{Symbol,1}:
 :val
 :adj



In [4]:

    
aa.val, aa.adj









    Out[4]:





(1.0,2.0)

Base.show(io::IO, x::Dual) = print(io, "$(x.val) + $(x.adj)ϵ")



In [5]:

    
aa









    Out[5]:





1.0 + 2.0ϵ

Making `Dual` work: conversion and promotion

To make Dual "just work" we need to define how to upconvert other numbers to it. We can do this by extending convert.

Base.convert(::Type{Dual}, x::Real) = Dual(x, 0) # first argument is not named, but is if type Type{Dual} Base.convert(::Type{Dual}, x::Dual) = x # corner case, since Dual <: Real

Notes:

the first argument is not named (nothing before the ::)
it is of Type{Dual}; Type{T} is a type of which the DataType T is the only member (that is, we're specifying tha that the call looked like convert(Dual, 4.) or similar
x can be any real number
we forward construction on to Dual; Real numbers just have adj = 0.



In [6]:

    
convert(Dual, 5)









    Out[6]:





5.0 + 0.0ϵ

We also need to define a promotion rule. Promotion rules stipulate what happens when we need to find a type "big enough" to hold two things (for addition, comparison, etc.)

Base.promote_rule(::Type{Float64}, ::Type{Dual}) = Dual Base.promote_rule(::Type{Int64}, ::Type{Dual}) = Dual Base.promote_rule(::Type{Float32}, ::Type{Dual}) = Dual

Arithmetic

With promotion and conversion defined, we only need to add arithmetic for Duals. Other arithmetic operations automatically promote as needed.

import Base: +, -, /, *, abs, <, ==, <=

+(x::Dual, y::Dual) = Dual(x.val + y.val, x.adj + y.adj) -(x::Dual) = Dual(-x.val, -x.adj) -(x::Dual, y::Dual) = x + -y *(x::Dual, y::Dual) = Dual(x.val * y.val, x.val * y.adj + x.adj + y.val) /(x::Dual, y::Dual) = Dual(x.val/y.val, x.adj/y.val - (x.val * y.adj)/y.val^2) abs(x::Dual) = x.val > 0 ? x : -x ==(x::Dual, y::Dual) = x.val == y.val && x.adj == y.adj <(x::Dual, y::Dual) = x.val == y.val ? x.adj < y.adj : x.val < y.val <=(x::Dual, y::Dual) = x < y || x == y



In [7]:

    
Dual(5, 2) - Dual(3, 1)









    Out[7]:





2.0 + 1.0ϵ



In [8]:

    
Dual(5, 2)/Dual(3, 1)









    Out[8]:





1.6666666666666667 + 0.11111111111111105ϵ

Now watch what we get for free:



In [9]:

    
5 + Dual(3, 2), Dual(3, 2) + 2









    Out[9]:





(8.0 + 2.0ϵ,5.0 + 2.0ϵ)



In [10]:

    
4.5/Dual(1, 1)









    Out[10]:





4.5 + -4.5ϵ



In [11]:

    
Array{Dual}(rand(5, 5))









    Out[11]:





5x5 Array{Dual,2}:
 0.08117297108120636 + 0.0ϵ  …  0.7761537892483081 + 0.0ϵ 
 0.05062865704382702 + 0.0ϵ     0.3574831527572737 + 0.0ϵ 
 0.7042205664707128 + 0.0ϵ      0.7852508128888431 + 0.0ϵ 
 0.7750474721966927 + 0.0ϵ      0.45371424589626574 + 0.0ϵ
 0.7305473566017255 + 0.0ϵ      0.7465167317384658 + 0.0ϵ



In [12]:

    
A = rand(5, 5) + Dual(0, 1) * Array{Dual}(rand(5, 5))









    Out[12]:





5x5 Array{Dual,2}:
 0.14552328892047406 + 1.5999828984728819ϵ  …  0.6566184836032978 + 1.28415435665137ϵ   
 0.9110217627851693 + 1.1730441782246503ϵ      0.5056012271322503 + 1.1224056065104175ϵ 
 0.8297085162113442 + 1.5845091401000588ϵ      0.886713630546704 + 1.3861618580595891ϵ  
 0.13998863884048807 + 1.210252817476958ϵ      0.21119856573946483 + 1.5117890146065958ϵ
 0.8559721429778611 + 1.1646014444060433ϵ      0.6604833718120571 + 1.9889443641540232ϵ



In [13]:

    
A * rand(5)









    Out[13]:





5-element Array{Dual,1}:
 0.3467906904838485 + 8.156917884377215ϵ 
 0.5993701684355719 + 7.633324292347406ϵ 
 0.6716694785307615 + 8.483065746839014ϵ 
 0.24579398447062234 + 7.725375280944693ϵ
 0.7679065688429951 + 8.526068040812186ϵ



In [14]:

    
minimum(A), maximum(A)









    Out[14]:





(0.03875850773520528 + 1.6630320021413ϵ,0.9633881400461819 + 1.1942231446682094ϵ)



In [15]:

    
Dual(1, 1) <= Dual(1, 0.8)









    Out[15]:





false



In [16]:

    
A \ rand(5)









    Out[16]:





5-element Array{Dual,1}:
  0.11507610919563258 + -118.12697765176911ϵ
  0.09480687542650598 + 956.749758164928ϵ   
  0.8384680698058897 + -939.7213129761792ϵ  
  0.45322362394265586 + -1038.554659367485ϵ 
 -0.43926866649295715 + 216.8425258556168ϵ

Automatic derivatives!

import Base.log log(x::Dual) = Dual(log(x.val), x.adj/x.val)



In [17]:

    
function foo(x, y)
    log(x + y) - y
end









    Out[17]:





foo (generic function with 1 method)



In [18]:

    
foo(Dual(1, 1), 2)









    Out[18]:





-0.9013877113318904 + 0.3333333333333333ϵ



In [19]:

    
foo(1, Dual(2, 1))









    Out[19]:





-0.9013877113318904 + -0.6666666666666667ϵ



In [ ]:

Project: Defining new types

Notes:

Making Dual work: conversion and promotion

Arithmetic

Automatic derivatives!

Making `Dual` work: conversion and promotion